Phase space factors for double-/? decay 
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A complete and improved calculation of phase space factors (PSF) for 2^/3/3 and Qvf5j3 decay is 
presented. The calculation makes use of exact Dirac wave functions with finite nuclear size and 
electron screening and includes life-times, single and summed electron spectra, and angular electron 
correlations. 
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I. INTRODUCTION 



Double-/? decay is a process in which a nucleus (A, Z) decays to a nucleus (A, Z ± 2) by emitting two electrons (or 
positrons) and, usually, other light particles 



(A, Z) (A, Z ± 2) + 2e T + anything. 



(1) 



Double-/3 decay can be classified in various modes according to the various types of particles emitted in the decay. 
For (3-/3-, the process 2v(3/3, Fig. QJi, 
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FIG. 1: Double-/3 decay mechanism for (a) two-neutrino, (b) neutrinoless and (c) neutrinoless decay with Majoron emission. 



(A, Z) -> (A, Z + 2) + 2e- + 2D 



(2) 



is allowed by the standard model and expected to occur with calculable probability. In recent years, the process Oz^/3/3, 
Fig. lb, 



(A, Z) -> (A,Z + 2) + 2c 



(3) 



has become of great interest, due to the discovery of neutrino oscillation [H-Q- The process is of utmost importance 
for obtaining the neutrino mass since its decay probability is proportional to the square of the average neutrino mass 
(m u )- A third process has been also considered, 0v(3f3M, Fig. [Tfc, 



(A, Z) (A, Z + 2) + 2er + M° 



(4) 
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in which a massless Nambu-Goldstone boson, called a Majoron, is emitted. However, most of the interest in this mode 
has disappeared in recent years and hence it will not be considered here. For /3 + /3 + decay, the corresponding modes 
2z//3/3, Oz//3/3, are 

(A, Z) -> (A, Z - 2) + 2e+ + 2v 
(A,Z) -> (A,Z-2) + 2e+. 

In this case, there are also the competing modes in which either one or two electrons are captured from the electron 
cloud, 2vfiEC, 2vECEC, Qv(3EC, QvECEC. 

For processes allowed by the standard model (2is/3f3, 2vj3EC ', 2vECEC) the half-life can be, to a good approxima- 
tion, factorized in the form 

= G 2 „\M 2l ,\ 2 , (6) 

where G% v is a phase space factor and M% v the nuclear matrix element. For processes not allowed by the standard 
model the half-life can be factorized as 

1 = Go„\M „\ 2 \f( mi ,U ei )\\ (7) 

where Gq v is a phase space factor, M$ v the nuclear matrix element and /(to,-, U e i) contains physics beyond the 
standard model through the masses and mixing matrix elements U e i of neutrino species. For both processes, two 
crucial ingredients are the phase space factors and the nuclear matrix elements. Recently, we have initiated a program 
for the evaluation of both quantities. For the nuclear matrix elements we have developed an approach based on the 
microscopic interacting boson model (IBM-2) and presented some results in Additional preliminary results have 
been presented in [a, Q and will be discussed in a forthcoming publication pj . In this article, we concentrate on phase 
space factors. 

A general theory of phase s pac e factors in DBD was developed years ago by Doi et al. [1, Q following previous 

work of Primakoff and Rosen jlOl | and Konopinski fill ] . It was reformulated by Tomoda (121] whose work we follow 
here. Tomoda also presented results in a selected number of nuclei. These results were obtained by approximating 
the electron wave functions at the nuclear radius and without inclusion of electron screening. In this article we take 
advantage of some recent developments in the numerical evaluation of Dirac wave functions and in the solution of the 
Thomas-Fermi equation to calculate more accurate phase space factors for double-/? decay in all nuclei of interest. 
Our results are of particular interest in heavy nuclei, aZ large, where relativistic and screening corrections play a 
major role. Studies similar to ours were done for single-/? decay in the 1970's [13] . In this article we report results for 
/3~j3~, which at the moment is the most promising decay mode. In a subsequent publication, we will present results 
for /3 + /3 + , (3 + EC, ECEC, which is very recently attracting some attention [14| . 
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II. ELECTRON WAVE FUNCTIONS 



The key ingredients for the evaluation of phase space factors in single- and double- j3 decay are the (scattering) 
electron wave functions. (For EC the bound wave functions.) The general theory of relativistic electrons can be found 
e.g., in the book of Rose fl5| |. We use, for (3 decay, positive energy Dirac central field wave functions, 

A K ,(r)=(ji\ r M V (8) 



where Xk are spherical spinors and g K (e, r) and f K (e, r) are radial functions, with energy e, depending on the relativistic 
quantum number k defined by n — (I — j)(2j + 1). Given an atomic potential V(r) the functions g K (e, r) and f K (e, r) 
satisfy the radial Dirac equations: 

dg K (e,r) k e - V + m e c 2 

dr r ch ^ 

df K {e,r) e-V-m e c 2 k 

— -i = z g K {e,r) + -f K (e,r). 

dr ch r 

The electron scattering wave function, denoted here by e s (e,r), where s is the projection of the spin, can then be 
expanded in terms of spherical waves as 

e s (e, r) = ef 1/2 (e, r) + ef 1/2 (e, r) + ef 3/2 (e, r) + ... (10) 
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where 



e Sl/2 ( e r)-f 9~i(e,r) X s 
e s ^*)-y h (z,r){v-a) X s 

e^(e,r)=( i9 ^f-^:^) (11) 

^3/2, s _ ( i#_ 2 (e,r)[3(r ■ p) - (r • a)(p ■ a)]x s 



er /2 (e,r) = 



«/2(e,r)[3(r • p)(p • a) - (r • a)] Xs 



The large and small components g K (e,r) and f K (e,r), respectively, with e — \/ (m e c 2 ) 2 + (pc) 2 of the radial wave 
functions are normalized so that they asymptotically oscillate with 




f*( e > r U ' W \ J i=^£ co8 (kr - q - r)ln(2kr) + 5 K ) 



(12) 



where 



n ch 

is the electron wave number, 77 = Ze 2 /hv is the Sommerfeld parameter and 8k is the phase shift. (For the neutrino 
wave functions appearing in the 2v decay mode the limit Z —> is taken, in which case the wave functions become 
the spherical Bessel functions.) 

The radial wave functions are evaluated by means of the subroutine package RADIAL [16], which implements a 
robust solution method that avoids the accumulation of truncation errors. This is done by solving the radial equations 
by using a piecewise exact power series expansion of the radial functions, which then are summed up to the prescribed 
accuracy so that truncation errors can be completely avoided. The input in the package is the potential V. This 
potential is primarily the Coulomb potential of the daughter nucleus with charge Zd, V(r) = —Zd(ahc)/r. As in the 
case of single-/? decay [13J we include nuclear size corrections and screening. 

The nuclear size corrections are taken into account by an uniform charge distribution in a sphere of radius R = 
TqA 1 / 3 with r = 1.2 fm, i.e. 



V(r) = 



Z d (ahc) r > R 

-Z d {ahc) , r<R 



(14) 



The introduction of finite nuclear size has also the advantage that the singularity at the origin in the solution of 
the Dirac equation is removed. (Other charge distributions, for example a Woods-Saxon distribution, can be used if 
needed.) 

The contribution of screening to the phase space factors was extensively investigated in single-/? decay pjlGjj]- The 
screening potential is of order Vs on Z^a 2 and thus gives a contribution of order a = 1/137 relative to the pure 
Coulomb potential Vc oc Z^a. We take it into account by using the Thomas-Fermi approximation. The Thomas-Fermi 
function (fi(x), solution of the Thomas- Fermi equation 



d 2 ^ _ ip 3 / 2 
dx 2 \/x 



(15) 



with x = r/b and 



where ao is the Bohr radius, is obtained by solving Eq. (| 15[) for a point charge Z c i with boundary conditions 

<P(0) = 1, 

I \ 2 (17) 
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This takes into account the fact that the final atom is a positive ion with charge +2. With the introduction of this 
function, the potential V(r) including screening becomes 



V(r) = (p(r) x 



z d (ahc) r > R 

-Z d (ahc) (^P 1 ) , r<R 



(18) 



This can be rewritten in terms of an effective charge Z e g = Zdf(r) where Z e g now depends on r. In order to solve 
Eq. (|15p . we use the Majorana method described in .1$] which is valid both for a neutral atom and a positive ion. 
The method requires only one quadrature and is thus amenable to a simple solution. It is particularly useful here, 
since we want to evaluate screening corrections in several nuclei. The Thomas-Fermi electron density is approximate, 
especially at the origin. However, the screening correction is only of order a relative to the Coulomb potential and 
the error on this small correction is therefore negligible. (A better method would be to do an atomic Hartree-Fock 
calculation and then fit the result to the expansion 



V(r) 



(-Z d (ahc)/r) ciiexp(-bix), 



(19) 



where x — r/b as in Eq. ([ToT) . However, it has been shown in single-/? decay that this method gives results comparable 
to the Thomas- Fermi approximation [18J, except in very light nuclei, Z < 8, which we do not discuss here.) We also 
do not consider radiative corrections to the phase space factors which are of order a 3 and thus negligible to the order 
of approximation we consider in this article. 

In order to show the improvement in our calculation as compared with the approximate solution used in the 
literature we show in Fig. [2] a comparison of the radial wave functions for 150 Nd decay, Zd = 62, at e = 2.0 MeV. 




FIG. 2: Electron radial wave functions g_i(e,r), f-i(e,r) (left panel) and /i(e,r), gi(e,r) (right panel) for Zd = 62, e = 2.0 
MeV and R — 6.38 fm (vertical line). The notations WF1, WF2, and WF3 correspond to leading finite size Coulomb, exact 
finite size Coulomb and exact finite size Coulomb with electron screening, respectively. 
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III. PHASE SPACE FACTORS IN DOUBLE-/? DECAY 
A. Two neutrino double-/} decay 

The 21/0/3 decay, Fig. QJi, is a second order process in the effective weak interaction. It can be calculated in a way 
analogous to single-/? decay. Neglecting the neutrino mass, considering only S-wave states and noting that with four 
leptons in the final state we can have angular momentum 0, 1 and, 2, we see that both + —> + and + — > 2 + decays 
can occur. We denote by Q/3/3 the Q-value of the decay, by the excitation energy in the intermediate nucleus, and 
by A the excitation energy with respect to the average of the initial and final ground states, 

A = -W + E N - Ej = i(Q M + 2m e c 2 ) + E N - El (20) 
The situation is illustrated in Fig. [3J 



W 



FIG. 3: Notation used in this article. The example is for Nd decay. 




1. 0+ -> 0+ 2u/3/3-decay 

The differential rate for 0+ 0^ 2i//3/3-decay is given by ( [Ml, ^ ) 

dW 2v = (a (0) + a (1) cos6» 12 ) w 2y dLo 1 de l d6 2 d{cos6 l2 ) (21) 

where ei and e 2 are the electron energies, lo\ and lo 2 the neutrino energies, 9\ 2 the angle between the two emitted 
electrons, and 

g^(Gcos 6> c ) 4 22 

w 2 y = ^Tfr w 1 w 2 (pic)(p 2 c)eie 2 - (22) 

The quantities and are a sum of the contributions of all the intermediate states and depend on the energy En 
of the intermediate state in the odd-odd nucleus and on the nuclear matrix elements M 2v . Introducing the short-hand 
notation 



(K 



N 



ei + uji + (E N ) - Ei e 2 +ui 2 + (E N ) - E r ' 

_ 1 1 [Z6) 

~ ei + cj 2 + (E N ) - Ei + e 2 +u>i + {E N )-Ei' 
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where (En) is a suitably chosen excitation energy in the odd-odd nucleus, one can write to a good approximation, 



,(°) = 



1 



A?\M 2 „\ 2 A 2 



((K N ) + (L N )f + ±((K N )-(L N )) 2 
((K N ) + (L N )) 2 -±((K N )-(L N )) 2 



(24) 



where M% v are the nuclear matrix elements and and are products of radial wave functions. Since Eq. (|24|) 
is an approximation to the exact expression, which is, however, of crucial importance for the separation of the decay 
probability into a phase space factor and a nuclear matrix element we have investigated the dependence of and 
on the energy (En)- Since (En) appears both in the denominator of Eq. (|2"4")l through (Kn) and (Ln) and in 
the numerator through A 2 — [Wo/2 + (En) ~ Ei] 2 , the dependence on (En) cancels almost completely, as already 
remarked years ago by Tomoda jl2j ]. and as it is shown by explicit calculation in the following paragraphs. 
The functions f[°^ and fff are defined as 



f (0) 

in 

f (i) 
in 



with 



ir^T + i/nr + i/ ir + i/i 
-2Rc[/- i - i /ri + .r 1 i/r 1 i. 

/n = /i(ei)/i(e 2 ), 
f-\=9-i(ei)h(e 2 ), 
/i _1 = /l(ei)5-i(£2). 



— 1 12 



(25) 



(26) 



The functions g~i(e) and /i(e) are obtained from the electron wave functions. We have used several ways to obtain 
<?_i(e) and fi(e) following an approach similar to that used in single-/? decay. We write 



w(r)g—i(e, r)r 2 dr, 
w(r)fi(e, r)r 2 dr. 



(27) 



In approximation (I) we use the weighing function w(r) = S(r — R)/r 2 in which case 



ff-i(e) = 9-i(e,R) 
f 1 (e)=f 1 (e,R) ' 



(I) 



(28) 



that is the electron wave functions are evaluated at the nuclear radius r — R. This is the simplest approximation 
and is commonly used in single-/? decay. We adopt it in this article. In approximation (II) we use the weighing 
function w(r) = 3/R 3 for r < R and w(r) = for r > R (an uniform distribution of radius R). This is not a good 
approximation, since the inner states cannot decay due to Pauli blocking and the decay occurs at the surface of the 
nucleus. Nevertheless, it is sometimes used. It essentially amounts to an evaluation of g_i(e) and /i(e) at a radius 
as one can show by explicitly evaluating 



3 f R 

9 - 1 ^ = R3j 9-i(t,r)r 2 dr 



Me) 
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(II) 



(29) 



ji(e,r)r 2 dr 



The third and most accurate approximation (III) is that in which the weighing function is the square of the wave 
function, R n i(r), of the nucleon undergoing the decay, 



9-1® 
Me) 



\R n i{r)\ g-i(e,r)r 2 dr 



\Rni(r)\ 2 fi(e,r)r 2 dr 



(III) 



(30) 
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By using harmonic oscillator wave functions and assuming that only one orbital is involved, the integrals in Eq. (|30|) 
can be easily evaluated. The approximation (III) essentially amounts to an evaluation of g_i(e) and /i(e) at a radius 
y/ (r 2 ) nl - For harmonic oscillator wave functions 



with 



one has 



b 2 = — ~ LOA^im 2 , (32) 
Muj 



(33) 



This approximation has the disadvantage that it must be done separately for each nucleus. Since in this paper we 
are seeking greater generality and do not wish to make a commitment to definite nucleonic orbitals, we make use of 
approximation (I). However, our computer program is written in such way as to allow the possibility of using Eq. (|30[) 
instead of Eq. (|28|) . Also in Sect. IIVI we study in a specific case, 110 Pd, where the transition is between lgg/ 2 and 
I57/2 orbitals, the error we make by using Eq. (|28p instead of Eq. (|30[) . 

All quantities of interest are obtained by integration of Eq. (|21[) . In the approximation described above, all quantities 
are separated into a phase space factor (independent of nuclear matrix elements) and the nuclear matrix elements. 
The two phase space factors are 

n A2 rQpp+irieC 2 rQpp+nrieC 2 —ei rQea— ei— 62 , , 
(0) _ AA / / / AO) 



3 In 2 

X 



F 2u — 0.0 / / / Al 



(34) 

((K N ) 2 + (L N ) 2 + (K N ) (L N )^j W2ydu>ide 2 dei, 



where W2 is determined as W2 



2A 2 
91n2 

x 



/j+m c c nQpp+m e c — ei pQpp—ei—ez 



2[(K N ) 2 + (L N ) 2 ) +5(K N ) (L N ) 



Ai) 
Ju 



w 2 vdoj 1 de2dei, 



(35) 



e i ~~ e 2 — W\- It has become customary to normalize these to the electron mass 



m e c . Also since the axial vector coupling constant §a is renormalized in nuclei it is convenient to separate it from 
the phase space factors and define quantities 



G 



n(0 
2// 



g\(m e c 2 ) 2 ' 



(36) 



These quantities are then in units of y . From these, we obtain: 
(i) The half-life 



(ii) The differential decay rate 



2u 
r l/2 



G 2 °Jg\ \m e c 2 M 2 J 



dW 2v dG 2 °J 
dei de\ 



(37) 



(38) 



where M 2v = g\ \m e c 2 M 2u \ . 

(iii) The summed energy spectrum of the two electrons 



dW- 



2// 



G?(ei + e 2 — 2m e c 2 



■M 2v 



dG 



(o) 

2// 



d(e\ + e 2 — 2m e c 2 ) 



(39) 
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(iv) The angular correlation between the two electrons 



<*( £ i) = (4°) 



dG^J/de 



We can evaluate the phase space factors G 2 ^ for any value (En)- The dependence of G%J on A = (Qfjp + 2m e c 2 )/2 + 
(En) — Ei is shown in Fig. |4Jd for the specific case of 110 Pd decay. We see that G^J depends mildly on A (< 1%) 
except very close to threshold (En) = 0, where the dependence is ~ 7%. A similar situation occurs for G$. We 
have done a calculation of G^J and G^J in the list of nuclei shown in Table Q] with A from Ref. [2(| or estimated by 
the systematics A — 1.12 A 1 / 2 MeV, which approximately represents the energy of the giant Gamow- Teller resonance 
in the intermediate odd-odd nucleus. The obtained G^J values are also shown in Fig. [5] where they are compared 
with previous calculations [2l|. These values of A are those estimated in the closure approximation and should be 
combined with the closure matrix elements 



M 2i ,~ ( .41) 
\9a) A f AGT' y ' 

where M£ = (0+\ Y, nn , T n T n , \0f) and M%J = (Q+\ Y, nn , r n T n ,a n ■ a n >\0f). Here A F is the closure energy for 0+ 
states in the odd-odd intermediate nucleus and it can be approximately taken as the energy of the isobaric analogue 
state. 

In recent years, it has been suggested, that in some nuclei, the lowest 1 + intermediate state dominates the decay. 
This is called the single state dominance hypothesis (SSD) 122 21 i . This situation is likely to occur in 96 Zr, 100 Mo, 
110 Pd, and 116 Cd, where protons occupy mostly the lgg/2 level and neutrons mostly the 1*77/2 level, and in 128 Te, 
where protons occupy mostly the 2d 5 / 2 level and neutrons mostly the 2c? 3 / 2 level, which are spin-orbit partners of 
each other. In the SSD model the energy (En) is that of the single state (En) — E x +. We have done a calculation 

of G^J and G^J for the nuclei mentioned above in the SSD case. This is also shown in Table U in columns 3 and 5. 
In this case, G^J and G%J should be combined with the matrix elements 

M GT = (0£||T+ ( T||l + )(l+||T+a||0+) 

2v ±(Qw + 2m e c 2 ) + E l +-E I ' 
Finally, using our program, one can evaluate the sum 

V- r W (0^||r+a||l+)(l+||r+ C 7||0+) 

V N hiQpp + 2m ^) + e n -Ej {6) 

if the individual GT matrix elements are known from a calculation, and a similar sum for Fermi matrix elements. In 
this case, there is no separation between 2vf3/3 phase space factors and nuclear matrix elements. 



Ml Mg T 



We also have available upon request for all nuclei in Table U the single electron spectra, summed energy spectra and 
angular correlations between the two outgoing electrons. As examples we show the cases of 136 Xe — > 136 Ba decay, 
Fig. [51 of very recent interest to EXO experiment [34[ and the case of 82 Se — > 82 Kr, Fig. [71 of interest to NEMO 
experiment [35|]. The use of our "exact" calculation makes a considerable difference as shown in Fig. [H For the 
SSD case there is a difference in the single electron spectra at small energies ei, as is shown in Fig. [S]for 110 Pd, and 
previously emphasized in Refs. (25l[26j. 



2. 0+ -¥ 0J 2v PP -decay 

The decay to the excited + state, 0^ (Fig. [3]), is also of interest. The phase space factor for this decay can be 
calculated using the formulas of the previous subsection, with Qpp replaced by 

Q^-E x (0+)=Q^(0t) (44) 



The results of this calculation are shown in Table [TT1 
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a) 



r 



e/) + =0.893(ll)MeV 



'Ag 



'Pel 



2^=2.01785(64) MeV 



\ 0^=2.8922(16) MeV 



'Cd 



b) 



1.55xl0- 19 F 



1.5 x itr 19 ■ 



^ 1.45 xKr 19 - 



1.4 x io- l9 y 



1.35 x l(T ly ■ 



1.3 x 10 



-19 . 



■ Closure energy A= 1.1 2A 1/2 =1 1.75 MeV 
• Single state dominance A= 1.893 MeV 



10 11 12 



A [MeV] 



FIG. 4: Panel a) Skeleton of the 110 Pd decay scheme. The ground state of the intermediate 110 Ag nucleus is 1 + leading to the 

lowest possible value for En to be E.+ = 0.0 MeV. 

1 i 

Panel b) Behaviour of the phase phase factor G^J as a function of A. The value obtained using single state dominance 
hypothesis, A = 1.893 MeV, is denoted by a red circle and the value obtained using A = 1.12 x 110 1/2 MeV= 11.75 MeV is 
denoted by a blue square. 

3. 0+ -)■ 2+ 2v/3/3- decay 

The half-life for 0+ 2+ 2z//3/3-decay is given by equations similar to those of sect. IIII A lpl . I20I l36j . The lepton 
phase space factor is now 



F, 



(0)0+ 



In 2 



AO) 



(45) 



x {(K N ) - (L N )) w 2 uduJide 2 dei, 



with Qppfif) = Qpp — E x (2i), (Fig. [3]), from which the life-time can be calculated 



M. 



(2+) 



(46) 
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Nucleus G^(10- 21 y- 1 ) G^ RRn (lQ- 21 y- 1 ) G^ICT 21 y" 1 ) G^J ^ n (lQ- 21 y" 1 ) Q^(MeV) i(MeV) issp(MeV) 





1000U. 




-nyou. 




Q.Z, I ZZul ^tU^t 1 


7 71 7^ 
(.tit 






/!£ 1 7 




9fi Q7 




9 O^QOfil (7\ a 


Q /1 1 1 h 




82 Se 


1596. 




-1075. 




9 QQ^I 9f / 901 ~1 

£i . is is U -L £i \ £i\J -L J 


08 h 




96 7 


£81 








^ ^^7(^)^0^ 
O.oOUo ( i^zoy ) 


1 n Q7 
iu.y / 


9 90^ 


100 Mo 


oouo. 


4134. 


-2263. 


97fi9 
-z / uz . 


*3 0' : i44n('1 7^ b 


11.20 


-L .yJOO 


no pd 


137.7 


146.9 


-79.56 


-84.45. 




11.75 


1.893 


116 Cd 


2764. 


3176. 


-1857. 


-2108. 


9 81 350f1 3~l d 


12.06 


1.875 


124 Sn 

Oil 


UOO.KJ 




"M9 7 




9 9&fiQ7M 


12.47 




128 T 


PI 9R88 

U.iUOO 


fl 9797 


-U. 1U1 i 


-U. 1UU1 


fl 8(^87^ 1 ^1 ^ e 

U.OUUO i ^lOl ) 


1 9 


1 fi&^i 


130 Te 


1529. 




-993.9 




9 ^9(\Q7(9^) d 


13.27 h 




136 Xe 






Q97 9 

C7i 1 .Z 




9 4^7S^('^7^ 


1 "3 flfi 




148 Nd 


324.8 




-195.5 






13.63 
















1 "3 79 
10. / z 




154 Sm 


9.591 




-4.816 




1.21503(125) 


13.90 




160 Gd 


193.8 




-114.2 




1.72969(126) 


14.17 




198p t 


15.36 




-8.499 




1.04717(311) 


15.76 




232 Th 


11.31 




-6.779 




0.84215(246) 


17.06 




238jj 


14.57 




-9.543 




1.14498(125) 


17.28 





TABLE I: Phase space factors G^J and G^/ obtained using screened exact finite size Coulomb wave functions. The Q-values 



-<(i) 



Ref. [30j|, 6 Ref. |3j, / Ref. [32} , 9 Ref. [33j) when available, 
1 or estimated by the systematics, A — X.YIA 1 ^ MeV, where 



are taken from experiment ( a Ref. [23], " Ref. [2j|, c Ref. [2 
or from tables of recommended values. A is taken from h Ref. @ 

A without tilde denotes the mass number. Phase space factors G^J ssu and G^J SSD correspond to values obtained using the 
SSD model, in which case the used j4ssd is listed in the last column. 



10000 



1^ 1000 



1 

o 



100 



Sis 10 



approximate 
this work 




Te 



40 60 80 100 120 140 160 180 200 220 240 

Mass number 



FIG. 5: Phase space factors G^J in units (10 21 y 1 ). The label "approximate" refers to the results obtained by the use of 
approximate electron wave functions. The figure is in semilogarithmic scale. 



The nuclear matrix elements can be written, in the closure approximation, as 

(47) 



.4 ;! 
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ei-m^pVIeV] e 1 +e 2 -2m e c 2 [MeV] e!-m e c 2 [MeV] 



FIG. 6: Single electron spectra (left panel), summed energy spectra (middle panel) and angular correlations between two 
outgoing electrons (right panel) for the 136 Xe — > 136 Ba 2i//3/3-deca,y. The scale in the left and middle panels should be multiplied 
by Miv when comparing with experiment. 




ei-m e c 2 [MeV] ei+e 2 -2m e c 2 [MeV] e!-m e c 2 [MeV] 



FIG. 7: Same as Fig.[6]for the 82 Se ^ 82 Kr 2^/3/3-decay. 

where 

m gt(2+) = (2 +u j2 TnTn , [ ffn ffn ,f) || +). (48) 

nn' 

Since this decay contains the term (if at) — (Ln), it is suppressed, due to cancellations, and it will not be considered 
further. Also, other models (SSD, no-closure) can be used, if needed. 

B. Neutrinoless double-/? decay 

The theory of fV/3/3 decay was first formulated by Furry [33] and further developed by Primakoff and Rosen , 
Molina and Pascual (36[, Doi et al. Q, and, Haxton and Stephenson [2(|. Here we follow mainly the formulation 
of Tomoda [l2[ • The phase space factors for Ov(3/3 decay are simpler than those of 2fj3/3 because of the absence of 
integration over the neutrino energies. Also, with two leptons in the final state and S-wave decay we can only form 
angular momentum 0, 1 and therefore the decay to 2 + is forbidden. 

1. 0+ -> 0^ 0v(3l3-decay 
The differential rate for the decay is given by 0, [l2| 

dW u = (a (0) + a (1) cos6»i 2 ) w ^eid(cos<9i 2 ) (49) 
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FIG. 8: Same as Fig. ©for the 150 Nd -> 150 Sm 2i//3/5-decay. The fi gure also shows the difference between our "exact" calculation 
and the previously used approximate calculation. 
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FIG. 9: Single electron spectra for the 110 Pd — > 110 Cd 2i//3/3-decay obtained using the two approximations discussed in the text, 
namely closure approximation and single state dominance hypothesis. 



where t\ and e 2 are the electron energies, 0i 2 the angle between the two emitted electrons, and 

g\{G cos B C ) A 



w 0v 



16tt 5 



(m e c ) (he )(pic)(p 2 c)eie 2 



(50) 



This decay is forbidden by the standard model and can occur only if the neutrino has mass and/or there are right- 
handed currents. In view of recent experiments on neutrino oscillations [H-0 it appears that neutrinos have a mass 
and we therefore consider the phase space factors for this case. The quantities a^ ' and a^ 1 ' in Eq. (|4"5)) can then be 
written as UM 



,(0 



f(i) 
ill 



to. 



I M 0w I 



(51) 



i = 0, 1, where Mq u is the nuclear matrix element and f[i\ f$ are the quantities given in Eq. (|25[) 



All quantities of interest are then given by integration of Eq. (|49|). Introducing 

2 pQi3/3+m c c 2 



F, 



(i) 



where e 2 is determined as e 2 = Qpp + m e c 



°" In 2 

2 e 1; and defining the quantities 



G, 



g\{AR?) 



(52) 



(53) 



where R = r^A 1 ^, r$ — 1.2 fm, is the nuclear radius, we can calculate: 
(i) The half-life 

2 



'1/2 



- r {0) n 4 

— U 0u 9 A 



{m v ) 



m, 



\M „r, 



(54) 
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Nucleus 


G^(10- 21 y- 1 ) 


G^IO" 21 y" 1 ) 


G«(10" 21 y" 1 ) 




£(0+)(MeV) 


Q^(0+)(MeV) 


48 Ca 


0.3627 




-0.1505 




2.99722(16) 


1.27504(253) 


76 Ge 


0.06978 




-0.02380 




1.122283(7) 


0.916757(167) 


96 Zr 


175.4 


185.3 


-103.8 


-109.2 


1.14813(7) 


2.20224(296) 


100 Mo 


60.55 


65.18 


-33.54 


35.89 


1.13032(10) 


1.90408(27) 


no pd 


0.004842 


0.004864 


-0.001371 


-0.001377 


1.47312(12) 


0.54773(76) 


116 Cd 


0.8727 


0.8878 


-0.3642 


-0.3701 


1.756864(24) 


1.056636(154) 


124 Sn 


0.01988 




-0.006408 




1.657283(22) 


0.629687(1552) 


iso Te 


0.07566 




-0.02705 




1.79352(11) 


0.73345(34) 


136 Xe 


0.3622 




-0.1451 




1.578990(23) 


0.878840(393) 


148 Nd 


0.009911 




-0.003339 




1.42446(4) 


0.50429(196) 


150 Nd 


4329. 




-2934. 




0.740382(22) 


2.630998(222) 


154 Sm 


0.01850 




-0.006583 




0.6806673(18) 


0.5343627(12518) 


160 Gd 


0.006318 




-0.002178 




1.279941(23) 


0.449749(1283) 




0.00004221 




-0.00001944 




0.69142(9) 


0.15073(255) 


238-jj 


0.0004635 




-0.0002289 




0.94146(8) 


0.20352(133) 



TABLE II: Phase space factors G%J and G%J for decay to the first excited + states, 0^, obtained using screened exact finite 
size Coulomb wave functions. Phase space factors G2^ SSD and G^'ssd correspond to values obtained using the SSD model. 



(ii) the single electron spectrum 

2fif(e 1 )w 0u (e 1 ) (55) 



dW 0u dG$> 

: Jvov — ; — -Nov 



dei dei 
where Af 0v = g\ \(m u )/m e \ 2 \M 0v \ 2 . 

(iii) and the angular correlation between the two electrons 

a(ei) . fj^i _ ^cg*. 

The factor (4J? 2 ) has been introduced in Eq. (|53|) to conform with standard notation [2l[, in which the nuclear matrix 
elements M 0l/ are given in dimensionless units, that is they are multiplied by R. The factor of 4, which is missing in 
Tomoda's definition but is necessary to make the calculation consistent with Boehm and Vogel, has been the cause of 
considerable confusion in the literature, as well as the value of r$ used in R = tqA 1 ^ 3 . Some authors use tq = 1.1 fm 
instead of tq — 1.2 fm. 

We have done a calculation of Gq*J and GqJ in the list of nuclei shown in Table IIII1 The obtained Gq*J values are 
also presented in Fig. [TU] where they are compared with previous calculations [2l[ . 

We also have available upon request the single electron spectra and angular correlation for all nuclei in Table Mil 
An example, 76 Ge decay, is shown in Fig. [TO 

2. 0+ -> Oj Ou/3/3-decay 

The decay to \ can also be calculated as in the previous subsection MI A 21 The results are shown in Table IIVI 



IV. EVALUATION OF THE ERROR 



The input parameters in the calculation of the phase space factors (PSF) are the Q-value, Qpp, and the nuclear 
radius, R. We take the Q value from experiments whenever possible and thus the error introduced in G is directly 
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Nucleus 


G'°'(10- 15 y" 1 ) 


G™(10- 15 y- 1 ) 


Q««(MeV) 


48 Ca 


24.81 


-23.09 


4.27226(404) 


76 Ge 


2.363 


-1.954 


2.03904(16) 


82 Se 


10.16 


-9.074 


2.99512(201) 


96 Zr 


20.58 


-18.67 


3.35037(289) 


100 Mo 


15.92 


-14.25 


3.03440(17) 


no pd 


4.815 


-4.017 


2.01785(64) 


116 Cd 


16.70 


-14.83 


2.81350(13) 


124 Sn 


9.040 


-7.765 


2.28697(153) 


i28 Te 


0.5878 


-0.3910 


0.86587(131) 


130 Te 


14.22 


-12.45 


2.52697(23) 


136 Xe 


14.58 


-12.73 


2.45783(37) 


148 Nd 


10.10 


-8.506 


1.92875(192) 


150 Nd 


63.03 


-57.76 


3.37138(20) 


154 Sm 


3.015 


-2.295 


1.21503(125) 


160 Gd 


9.559 


-7.932 


1.72969(126) 


198p t 


7.556 


-5.868 


1.04717(311) 


232 Th 


13.93 


-10.95 


0.84215(246) 


238-y 


33.61 


-28.13 


1.14498(125) 


TABLE III: Phase space 


factors GqJ and G^J obtained using 


screened exact finite size 


Coulomb wave functions. 




80 100 120 140 160 180 200 220 240 

Mass Number 

FIG. 10: Phase space factors Gq° in units (10~ 15 y _1 ). The label "approximate" refers to the results obtained by the use of 
approximate electron wave functions. The figure is in semilogarithmic scale. 



related to the experimental error. For example, recently the Q -va lue for 110 Pd decay has been measured with high 
accuracy [2f| . Table |V] shows the improvement in the error in G^J and G^J due to the better accuracy obtained by 
measurement compared to the Q- value determined from mass values. 



The nuclear radius enters in the calculation in various ways, the most important of which is the evaluation of the 
quantities <?_i(e) and /i(e). We evaluate the error here by comparing approximation (I) with (III) in a specific case, 
110 Pd, where the transition is I59/2 — I57/27 obtaining an estimate of the error of 3%. For 0V decay the radius R 
enters also in the definition of Gq v . This is, however, an input parameter which does not depend on the method of 
calculation. We have used R = tqA 1 / 3 with ro = 1.2 fm. We can estimate the error introduced by this choice by the 
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m e c z [MeV] 



FIG. 11: Single electron spectra (left panel), and angular correlations between the two outgoing electrons (right panel) for the 
76 Ge — > 76 Se Oz^/3/3-decay. The scale of the left panel should be multiplied by A/oi/ for a realistic estimate. 



Nucleus 




G«(io- 15 y- 1 ) 


£(0+)(MeV) 


Q w (0+)(MeV) 


48 Ca 


0.2989 


-0.2080 


2.99722(16) 


1.27504(253) 


76 Ge 


0.1776 


-0.09855 


1.122283(7) 


0.916757(167) 


96 Zr 


4.566 


-3.760 


1.14813(7) 


2.20224(296) 


100 Mo 


3.162 


-2.493 


1.13032(10) 


1.90408(27) 


uo pd 


0.08844 


-0.02958 


1.47312(12) 


0.54773(76) 


116 Cd 


0.7163 


-0.4075 


1.756864(24) 


1.056636(154) 


124 Sn 


0.1709 


-0.06237 


1.657283(22) 


0.629687(1552) 


i30 Te 


0.3086 


-0.1271 


1.79352(11) 


0.73345(34) 


136 Xe 


0.6127 


-0.2924 


1.578990(23) 


0.878840(393) 


148 Nd 


0.2010 


-0.05354 


1.42446(4) 


0.50429(196) 


150 Nd 


27.27 


-23.26 


0.740382(22) 


2.630998(222) 


154 Sm 


0.2806 


-0.07744 


0.6806673(18) 


0.5343627(12518) 


160 Gd 


0.2063 


-0.04650 


1.279941(23) 


0.449749(1283) 




0.2622 


-0.0.1065 


0.69142(9) 


0.15073(255) 


238-jj 


0.7534 


-0.03918 


0.94146(8) 


0.20352(133) 



TABLE IV: Same as Table UTTl but for the decay to the first excited 0+ state, 0% . 



same method used in the phase space factors for single-/? decay [l3j], that is by adjusting r for each nucleus, A, Z, 
using 

\rlA^ = (r*) exp , (57) 

where (r 2 ) exp is obtained from electron scattering and/or muonic x-rays. The largest difference between R t h and R eX p 
is found to be ~ 4%. This leads to an error estimate of 0.5% for 2v. For Qv we obtain an estimate of error of 7%. 

In addition, we have an error coming from screening and most importantly from the value of (En)- We estimate 
the screening error to be 10% of the Thomas-Fermi contribution, known to overestimate the electron density at the 

nucleus. This gives an error in 

G oJ> G 2v of 0.1%. The estimate of the error introduced by the choice of (En) is 
model dependent. If we vary A from the value 1.12A 1 / 2 MeV to the SSD value (~ 2 MeV) we obtain for 110 Pd decay 
an error of 7%, as shown in Fig. \^p. If, however, we stay within a specific model, closure or SSD, the error estimate 
is much smaller. In particular for the SSD model the error is only arising from the value of Qec and Qpp shown 
in Fig. [4ji. The estimate therefore depends on the nucleus considered. For 110 Pd, the SSD model appears to be a 
good approximation and using it we obtain an estimate of the error of 0.05%. For the closure approximation the 

dependence of G^j on A is very mild (< 1%) except very close to the threshold, (En) = 0, as shown in Fig. |Ud. The 
situation is summarized in Table IVII 
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Qpp keV 






2004.00(1133)° 
2017.85(64) b 


1.386(67) x 10~ 19 
1.469(05) x 10~ 19 


4.707(86) x 10" 15 
4.815(06) x 10~ 15 



TABLE V: The uncertainty on PSF due to the uncertainty of the Q value. a From Ref. [3£| and b from Ref. [2!|. 



2v Q- value 10 x SQ/Q 

Radius 0.5% 
Screening 0.10% 

(En) model dependent 



0v Q- value 3 x SQ/Q 

Radius 7% 

Screening 0.10% 

(E N ) 



TABLE VI: The estimate of uncertainties introduced to phase space factors G 2 J and GqJ due to different input parameters. 



V. USE OF PHASE SPACE FACTOR 



The main use of phase space factors (PSF) is in connection with a calculation of the nuclear matrix elements to 
predict life-times for the decay. Here an important point is that the nuclear matrix elements are defined in a way 
consistent with the phase space factors. For example, we have defined the phase space factors for Oi//3/3 with a factor 
of 4 in Eq. (|53l) . This factor is not included in Tomoda's definition [l2| but it is in the book of Boehm and Vogel 
[2lT |. The nuclear matrix elements consistent with this factor are, for GT, those of J2 n n' T n T n"^n ■ ov, not those of 
(1/2) J2 n n' T nTn'&n ■ 0V. We will present results of our predictions where phase space factors are combined with the 
IBM-2 nuclear matrix elements in a forthcoming publication [3{|. Here we use the calculation of PSF to extract the 
1v matrix elements from experiments where the life-time of 2z//3/3 decay has been measured. The quantity we extract 
is the dimensionless quantity g^|(m e c 2 )M 2l/ | 2 = |M|*| 2 (also called Ni v in Sect IIII ATT) . The extraction of |Mf^| is 
possible in two cases: (1) the closure approximation (CA) and (2) the single state dominance (SSD) hypothesis. If 
neither of these two approximations is valid, then the quantities G21, and M<i v cannot be separated as discussed after 
Eq. (|43|) . The results obtained with the assumption of CA with A = 1A2A 1 / 2 MeV and under the assumption of SSD 
for 96 Zr, 100 Mo, 116 Cd, and 128 Te are shown in Table IVTI1 and in Fig. [T^J We note that all effective matrix elements in 
Table IVTT1 vary between a minimum of ~ 0.02 ( 136 Xe) and a maximum of ~ 0.2 ( 100 Mo and 238 U), with the majority 
being ~ 0.05. 



The effective matrix elements |M^| cxp can, in principle, be obtained from measurements of GT^ strengths (and 
F ± strengths), through the formula 



,<.-/ _ V (0^||T+a||l+)(l+||r+q||0+; 



N 

|^2^|° XP = 9a \( m e c ' 2 )M2i7'\i an d similar formulas for the Fermi matrix elements. However, in experiments, only the 
magnitude of the individual GT matrix elements can be measured, not its sign. Furthermore, it must be decided to 
what N to stop the evaluation of the sum, and what value to use for g^. Therefore, theoretical models must be used 
to obtain iMf^l for GT ± strengths. A recent example is 150 Nd decay [42| (one should note that in this paper the 
denominator in the definition of 1v Gamow- Teller nuclear matrix element is different from Eq. (UHl by 2m e c 2 due 
to the use of atomic masses in the calculation of Qpp and — Ej), where M^J has been extracted under (i) the 
assumption that only the 1 + state at £ I ( 150 Pm)= 0.11 MeV contributes to the decay and (ii) that all states up to 
£ I ( 150 Pm)< 3.0 MeV contribute. The result is (i) M 2 ^, T (MeV- 1 )= 0.028 ±0.006 and (ii) M£ T (MeV" 1 )= 0.13 ±0.02. 
Multiplying by (m e c 2 ) = 0.511 MeV and g\ = 1.273 2 0, one obtains (i) \M^\ C ^ D = 0.023 ± 0.005 and (ii) 
I M cff| ox P _ 108 ± 0.017. These two estimates bracket our extracted value 0.058 ± 0.004. This "experimental" way 
of extraction also assumes that the factorization of r 2 ^ 2 to G21/ and M-i v is valid. 
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Nucleus 


G'^IO-V" 1 ) 




^(lO 18 y) exp a 


|M 2 f | 


|M 2 f |ssd 


48 Ca 


15550. 




44 +l 


0.038 ± 0.003 




76 Ge 


48.17 




1500 ± 100 


0.118 ±0.005 




82 Se 


1596. 




92 ±7 


0.083 ± 0.004 




96 Zr 


6816. 


7825. 


23 ±2 


0.080 ± 0.004 


0.075 ± 0.004 


100 Mo 


3308. 


4134. 


7.1 ± 0.4 


0.206 ± 0.007 


0.185 ±0.006 


100 Mo- 100 Ru(0+) 


60.55 


65.18 


590l«° 


0.167 ±0.011 


0.161 ±0.010 


116 Cd 


2764. 


3176. 


28 ±2 


0.114 ±0.005 


0.106 ±0.004 


128 Te 


0.2688 


0.2727 


1900000 ± 400000 


0.044 ± 0.006 


0.044 ± 0.006 


128 Te 


0.2688 


0.2727 


3500000 ± 2000000 6 


0.033 ±0.017 


0.032 ±0.017 


130 Te 


1529. 






0.031 ± 0.004 




" 6 Xe 


1433. 




2110 ± 250 c 


0.0182 ±0.0017 




150 Nd 


36430. 




8.2 ±0.9 


0.058 ± 0.004 




150 Nd- 150 Sm(0+) 


4329. 




1331JS 


0.042 ± 0.006 




238 jj 


14.57 




2000 ± 600 


0.19 ±0.04 





TABLE VII: Experimental 2^/3/3 half-lives and the corresponding effective nuclear matrix elements |M|^|. For the case 128 Te 



two experimental half-lives are listed, upper from evaluation of Barabash a 
The value for 136 Xe is from a new measurement and is taken from 



and the lower from the comment of Pritychenko 
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FIG. 12: Effective nuclear matrix elements |M|* | extracted from the experimental 2uf3/3 half-lives as a function of mass number. 



Our calculation of G 2 J SSD allows one to test the SSD assumption for 100 Mo, 116 Cd, and 128 Te, where the matrix 
elements (even-even — > odd-odd) + — > 1± and (odd-odd — > even-even) 1^ — > + arc known from single (3 decay 
experiments. The extracted values of l-Ml^r** using the single j3 decay (or EC) matrix elements and gA = 1.273 are 
|M|5| CXP = 0.174 ± 0.075, 0.148 ± 0.023, 0.0152 ± 0.0003 for 100 Mo, 116 Cd, and 128 Te, respectively. These values, as 
well as that of 150 Nd discussed above, are given in Table IVllTl and compared with the ones obtained from experimental 
double j3 decay half-lives. The SSD model appears to give a rather good agreement for 100 Mo and 116 Cd, but is off 
by a factor of 2 in defect in 128 Te and 150 Nd. The situation has been also analyzed in detail from different point of 
view in Ref. 0. 
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Nucleus \Mg\ |AggJssp \M^\ a s f D 

100 Mo 0.206 ± 0.007 0.185 ±0.006 0.174 ±0.075 

100 Mo- 100 Ru(0^) 0.167 ±0.011 0.161 ±0.010 0.104 ±0.045 

116 Cd 0.114 ± 0.005 0.106 ±0.004 0.148 ±0.023 

128 Te 0.044 ± 0.006 0.044 ± 0.006 0.0152 ± 0.0003 

128 Te 0.033 ± 0.017 0.032 ± 0.017 0.0152 ± 0.0003 

150 Nd 0.058 ± 0.004 0.023 ± 0.005 



TABLE VIII: Effective nuclear matrix elements |Mf„ I, |AfsL |ssd obtained from experimental 2v/3f3 half-lives compared with 
lAf^lssD' the effective nuclear matrix elements obtained from single /3 decay experiments ( 100 Mo, 116 Cd, 128 Te), or from GT ± 
strength measurements ( 150 Nd). For the case 128 Te the two values listed are explained in the caption of Table fvTTI 



VI. CONCLUSIONS 



In this article, we have reported a complete and improved calculation of phase space factors for 2isf3~ f3~ and 
0vj3~ j3~ decay, including half-lives, single electron spectra, summed electron spectra, and electron angular correla- 
tions, to be used in connection with the calculation of nuclear matrix elements. Apart from their completeness and 
consistency of notation, we have improved the calculation by using exact Dirac wave function with finite nuclear 
size and electron screening. The program for calculation of phase space factors has been set up in such a way that 
additional improvements may be included if needed (P-wavc contribution, finite extent of nuclear surface, etc.) and 
that it can be used in connection with the closure approximation, the single state dominance hypothesis and the 
calculation with sum over individual states. In a subsequent publication we are planning to present complete and 
improved calculations for 2i//3 + /3 + and OV/3 + /3 + decay, as well as of the competing processes 2vf3 + EC, 2vECEC and 
Ovfi+EC, OvECEC. 
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